Generalized two-body self -consistent theory of random linear dielectric composites: 
an effective-medium approach to clustering in highly-disordered media 
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Effects of two-body dipolar interactions on the effective permittivity/conductivity of a binary, symmetric, ran- 
dom dielectric composite are investigated in a self-consistent framework. By arbitrarily splitting the singularity 
of the Green tensor of the electric field, we introduce an additional degree of freedom into the problem, in the 
form of an unknown "inner" depolarization constant. Two coupled self-consistent equations determine the latter 
and the permittivity in terms of the dielectric contrast and the volume fractions. One of them generalizes the 
usual Coherent Potential condition to many-body interactions between single-phase clusters of polarizable mat- 
ter elements, while the other one determines the effective medium in which clusters are embedded. The latter is 
in general different from the overall permittivity. The proposed approach allows for many-body corrections to 
the Bruggeman-Landauer (BL) scheme to be handled in a multiple-scattering framework. Four parameters are 
used to adjust the degree of self-consistency and to characterize clusters in a schematic geometrical way. Given 
these parameters, the resulting theory is "exact" to second order in the volume fractions. For suitable parameter 
values, reasonable to excellent agreement is found between theory and simulations of random-resistor networks 
and pixelwise-disordered arrays in two and tree dimensions, over the whole range of volume fractions. Com- 
parisons with simulation data are made using an "effective" scalar depolarization constant that constitutes very 
sensitive indicator of deviations from the BL theory. 

PACS numbers: PACS numbers: 05.60.Cd, 72.80.Tm, 78.20.Bh 



-a 

G 

O 
O 



> 

in 
oo 

o 
o 

(N 



X 



I. INTRODUCTION 



The problem of the effective transport properties of disor- 
dered media has a long history yj, |2J], and still continues to 
attract wide attention owing to its intrinsic theoretical interest 
|0-@|, and to its importance to various domains of engineer- 
ing sciences [ 6] . Our focus here is on the effective permittivity 
of binary symmetric random media H with quenched disor- 
der, that undergo percolative behavior [0 0]. Investigations 
have mainly been carried out on prototypical models such as 
Random Resistor Networks (RRN) [9-11], random arrays of 
polarizable point elements 01211 and checkerboards IIOl . The 
advent of full-field numerical methods of computation, such 
as Fast Fourier Transform (FFT) calculations on pixel arrays 
filEl (see fill for an alternative use of FFTs in this con- 
text), or finite-elements methods, make it now possible to ad- 
dress in detail random checkerboards and various other types 
of heterostructures 01711. 



From a theoretical standpoint, accounting for presence of 
a percolation threshold in a systematic theory of interac- 
tions between heterogeneities is extremely difficult, since the 
highly-disordered character of the percolation regime involves 
many -body positional correlation functions to all orders |00]. 
There, perturbative methods are inapplicable unless some 
amount of self-consistency is injected in suitable approxima- 
tions. The simplest successful [18] self-consistent (s.c.) ap- 
proach to the effective permittivity problem in percolating 
media is the well-known Bruggeman-Landauer (BL) theory 
filEU], whose theoretical status is well-established with re- 
gard to perturbative approaches O20l l2lll . It can be applied to 



RRNs and to continuum systems. In the BL theory, as well 
as in more elaborate s.c. treatments 1I20I I22I . the threshold is 
fixed, which is a nuisance since it varies in practice with mi- 
crostructure [|5|. 

The present work proposes a parametric theory to account 
for corrections to the BL theory in the high-contrast, highly- 
disordered regime. Simple ways of doing so mostly rely on 
varying the shape of the inclusions by modifying their de- 
polarization coefficients J23|[. The related formulations of 
McLachlan and co-workers [24] empirically introduce tunable 
critical exponents and threshold in the BL effective-medium 
formula. In contrast, our purpose is to incorporate corrections 
to the BL theory by considering exact pairwise interaction 
terms between polarizable pointlike matter elements, which 
will be done in the generic continuum framework of Miller's 
cell-material model ||25[ l26ll . In this connection, it should 
be mentioned that Shen and Sheng previously accounted for 
nearest-neighbor and next-nearest-neighbor interactions - a 
special case of pairwise terms - in a BL-like framework, by 
using two-dimensional split inclusions [27], which resulted 
in a theory with a double percolation threshold 02811 . This 
double-threshold effect is nowadays actively discussed l29ll . 
notably in the context of random checkerboards [30]. 

Our starting point is the modified BL equation for the ef- 
fective permittivity e e 



4e + (l-4)f 



0. 



(1.1) 



The brackets denote an average over material cells of vari- 
able permittivity e, and t e is an unknown effective depolariza- 
tion coefficient to be determined. In the original BL theory, 
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£ e = l/d where d is the space dimension, which is equal to 
the percolation threshold. Eq. ( Il.lt arises in a variety of con- 
texts. In particular, the Wu-McLachlan formula reduces to 
the above one when its exponents are set to one II24I1 . Rather, 
our line of thought will be that i e is a dimensionless func- 
tion of permittivity ratios and of the volume fraction of phases 
MM- 

The phenomenology of £ e is established in Sec. |IT]by iden- 
tifying the solution of Eq. (II . lb with results of simulations 
of binary RRN and pixelwise-disordered systems in two and 
three dimensions. In Sec. [Till we develop a formalism rooted 
in multiple-scattering theory, in which a free, inner, depo- 
larization constant £, analogous to the above effective one, 
can be introduced without making approximations. Indepen- 
dently, microstructure-related features are approximately ac- 
counted for by means of two parameters, which provides a 
rough characterization of single-phase clusters. With this sim- 
ple device, we make an bypass the need to consider explicitly 
correlation functions. Section ITVl formulates the two s.c. con- 
ditions required to determine £ and the effective permittivity. 
One of these modifies the usual coherent-potential condition 
of vanishing self-energy, into one that distinguishes between 
the local and non-local parts of the two-body term in the self- 
energy. It serves to adjust, by means of two supplementary 
parameters, the "amount of locality and non-locality" that en- 
ter the s.c. condition. This peculiar treatment is justified by 
comparing the resulting four-parameter theory with the simu- 
lations of Sec. [TO We conclude in SeclVI 



II. PHENOMENOLOGICAL BEHAVIOR OF l e 

A. Preliminary remarks 

We consider a <i-dimensional binary medium, with phases 
of permittivities E\ and £2 in respective proportions 1 — / and 
/. For definiteness, we choose phase 2 as that of high permit- 
tivity. Eq. d 1 - 1 b then reduces to the second-order polynomial 
equation for the effective permittivity e e 



(1-/)- 



E\ — e e 



+ f- 



£2 — £e 



' lex + (1 - £e)e e ' " ie 2 + (1 - te)e' e °' (2 ' !) 

In the dilute limit, an exact well-known HH] result (due to 
Maxwell in three dimensions) is that: 



e e / £l =1 + /- 



d{e 2 - si) / , , 



£2 + (d - 

On the other hand, expanding ( 12. It provides: 

£2 - El 



0(f). (2.2) 



e e /ei = 1 + / 



fe 2 + (l-4)ei 



0(f 



(23) 



where £ e stands for £ e (f = 0). The 0(f) term arises from the 
exact polarizability factor that individually characterizes the 
impurities (one-body term). Thus, £ e (f = 0) = l/d and the 
O(f) correction to this value is due to the 0(f 2 ) term in e e , 
which accounts for pairwise (two-body) interactions. Similar 
considerations hold near f — 1 where £ e (f = 1) = l/d. 



Next, percolative behavior takes place when E\ <C e e <C 
£ 2 . Letting e e /s\ — > 00 and £2/^1 00 at the percolation 
threshold / = f c in ( 12. Il l provides the following equation for 

fcM: 



Ufc) = fc 



(2.4) 



In the critical region, differences are expected between two- 
and three-dimensional cases. For infinite contrast, the mod- 
ified BL model ( 12.11) reduces to e e = £i(l — f/£ e ) for 
/ < / c , and e e = e 2 (f - £ e )/(l - 4) for / > f c . On the 
other hand, the well-known critical behavior of binary media 
in the infinite-contrast limit is [5] e e oc £i(l — /// c ) _s for 
f < f c and £ e oc £ 2 (/// c - 1)' for / > / c , with d-dependent 
critical exponents s and t. Comparing these equations pro- 
vides 

4(/) =i / + a-f c (l - f/f c ) s (f < f c ), (2.5a) 

~ / - o+(l - f c )(f/fc - If (f > / c ).(2.5b) 

where a ± > are coefficients of order one. For d — 2, ex- 
ponents are such s = t~1.3>l, so that derivatives at f c 
are 



^e(/c ± ) = l (d=2), 



(2.6) 



and £ e (f) is smooth. For d = 3 instead, exponents are now 
,s < 1 and t > 1 and 



-OO, 



an 



1 



(d = 3), 



(2.7) 



so that £ e (f) must display a cusp at f c . 

These considerations must be modified when applied to 
a mean-field theory such as the one developed hereafter, in 
which the critical exponents would have typical effective- 
medium values s = t = 1 irrespective of d. Then, 



aft) 



) = l + a+(l-/- 1 ). 



(2.8a) 
(2.8b) 



Depending on the values of a ± and / c , £' e (ft) can men ( a 
priori) be of either sign. 

To close, we point out that the double-threshold effect for 
two-dimensional checkerboards (see Introduction), translates 
in terms of £ e into the property £ e (f) — / for / G (p c ,l—p c ), 
where p c is the two-dimensional site-percolation threshold 
& 



B. Simulation data 

The BL theory compares well to binary RRNs, provided 
that an identification is made between the phase volume frac- 
tion and the volume fraction of bonds, because both are sym- 
metric with respect to phase interchange 112 U 13411 . In particu- 
lar, for d = 2 the percolation threshold coincide, which leads 
to an overall reasonable match. Moreover, the effective con- 
ductivity (permittivity) in the BL theory coincides with that 
of RRNs up to the third order (included) of perturbations in 
powers of the scaled contrast 5e/ (e) 12111 . and it is exact to 
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FIG. 1 . (Color online) Permittivity (conductivity) data analyzed in terms of l e (/) for two- and three-dimensional models, and contrasts £2 /ei 
as indicated in the plots: (a) and (b), random-resistor networks of size L — 512 (resp., L = 64) for d = 2 (resp., d = 3); (c) and (d), 
pixelwise-disordered media of same sizes (see text). Dashed: 1 jd. Quasi-vertical grey line: function £ e = /. Solid lines between data points 
are guides to the eyes. 



one-body order (included) in the sense of Sec. lIIICl However, 
strong discrepancies with RRNs arise for d > 3, primarily be- 
cause the percolation threshold of the BL theory (f c = l/d) 
overestimates that of the lattice system. 

Discrepancies can be analyzed by expressing £ e (f) as a 
function of e e by means of ( 12.11 ). as IB5I1 



£e«e) 



(e e - £l)(£2 - Ee)' 



(2.9) 



and by using for e e effective-permittivity data obtained from 
models. We consider three numerical models. 

To begin with, Fig. 02a) and (b) show an analysis of the ef- 
fective conductivity of simulations of bond-disordered RRNs, 
carried out for the present purpose. The function £ e (/) is com- 
puted from d2.9t for various contrasts. Statistical averages of 
the overall conductivity of a sufficiently large number of sam- 
ples have been carried out to reduce standard deviations to the 
typical error bar values represented in Fig. 03 a) for / < 1/2. 
Although not represented, they are of same order of magni- 
tude in the domain / > / c , and in Figs. |Tfb), (c) and (d). 
Under-sampling of the dilute configurations makes the statis- 
tical errors larger in the limits / — >• 0,1. System sizes are 
large enough to make finite-size effects negligible to our pur- 
pose. According to (12. 4K f c is the volume fraction / at the 
crossing point between the infinite-contrast plots and the grey 
line that represents £ e = f. It is seen that f c (d = 2) = 1/2, 
and f c (d = 3) ~ 0.26, close to the expected theoretical bond- 
percolation value ~ 0.249. The dilute-limit value £ e = l/d is 
represented by the dashed horizontal line, and is approached 
with negative slopes in all cases. Note that in three dimensions 
i e slightly exceeds 1/3 for / > 0.55. In this high-conductivity 
region, the closeness to Bruggeman's theory is remarkable. 

We also carried out noiseless calculations using 
Bernasconi's Real-Space Renormalization-Group model 
for RRNs IB6I1 . Up to some limitations of the approach, and 
to a lower percolation threshold in three dimensions, results 
are good qualitative agreement with Figs. [Tfa) and (b) [? ]. 

Figs. [Tfc) and (d) display further results of simulations 
on pixelwise-disordered arrays (PDA) 01 51 . solved using the 
FFT method developed for elastic composites by Michel et al. 
OI4I1 . Its adaptation to linear dielectric media is straightfor- 



ward. The formulation employed is the original one, that uses 
the continuum Green function (other types of implementation 
were considered in IU5I1 ). Each sample is a regular square 
or cubic array of L d pixels (voxels, in three dimensions), of 
permittivity chosen at random according to the binary prob- 
ability density. Conceptually closer to a random array with 
substitutional disorder, than to a bond network or a random 
checkedboard (the fields are not resolved within the pixels or 
voxels), this system nonetheless features in two dimensions 
the same percolation threshold f c = 1/2 as a bond-disordered 
RRN; for d = 3, f c ~ 0.315, a value reminiscent of site per- 
colation (f c ~ 0.312), close to the Bruggeman value. For 
d = 2 its £ e (f) function resembles that of RRNs, with a 
sharper variation at threshold. The situation changes in three 
dimensions, especially in the high-permittivity phase where 
deviations from the dilute limit markedly differ from that in 
RRNs: the approach of / = 1 has a positive slope; the con- 
cavity at small / is opposite; moreover, in the infinite-contrast 
limit, £ e (f) develops a valley just after the percolation thresh- 
old. The cusp at / ~ 0.36 marks out a transition from the 
critical region where Eq. (12.7) (right) applies, to behavior of 
the effective-medium type ( |2.8bt . The critical region / < f c 
where Eq. (12.7b (left) would lead to a cusp is not observed. 
This suggests that the critical behavior of this system, which 
has no contact interactions, might be different from the RRN 
one, with either s ~ 0, or s = 1 with a~ <C 1. This point 
-not crucial to our purpose- has not been investigated further, 
due to some difficulties in achieving high-quality numerical 
convergence for / < f c in the infinite-contrast limit (this is 
the reason why contrast is limited to 10 4 ). 



The rich typology of £ e (f) behaviors, even in the dilute re- 
gion, illustrates the dramatic influence of microstructural fea- 
tures on the overall response. This makes £ e {.f) an interesting 
means of analysis, since plots of £ e provide nontrivial infor- 
mation over the whole range of concentrations (trying to use 
it at low contrast on noisy data may however lead to an incon- 
clusive outcome 113 ill ). 
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III. THEORY 



Introducing moreover a "screened" electric field 



The theoretical framework we adopt to account for the 
above observations hea vily relies on the well-known multiple- 
scattering formalism ll37l HU, widely used in solid-state 
physics 39], and repeatedly employed to study dielectric me- 
dia Because some additions to the classical setup 
are needed, the following sections review it briefly with em- 
phasis on our modifications. A number of the equations also 
arise in the context of dilute alloys, a related problem, where 
Bruggeman's EMA is known as the Coherent-Potential Ap- 
proximation (CPA) i43ll : see l44ll for a recent review. 



A. Split-up of the Green function 

We consider a o?-dimensional dielectric medium with per- 
mittivity e(r) fluctuating from cell to cell. The latter are spher- 
ical on average, of infinitesimal typical radius a —> and vol- 
ume v = a S d /d, where S d = 2tt^ 2 /T(d/2) is the area of 
the d-dimensional sphere of unit radius. We use the notations 
I for the d x d identity matrix, and r = r/r for the unit po- 
sition vector. Following the usual treatment, we introduce an 
arbitrary background permittivity Eb, and re-write the equilib- 
rium equation V.(eE) = 0, where E is the electric field, as 
e&V.E = — V.fcE, where Se — e — £b is the dielectric con- 
trast. The problem can then be cast in the form of an integral 
equation for E P45I1 : 

E(r) = E (r) + / cfV' G 1/d (r - r')^E(r'), (3.1) 
J £b 

where Eo is the applied field, and G\/d is the d-dimensional 
dipolar Green function (6 is the Heaviside function, used to 
implement a principal value prescription at the origin) [46] 



I 



1 1 



G 1/d (r) = --5(r) - lim 6{r - rj) — -3(1 - dvv). (3.2) 



7^0 



In operator notation where, e.g., e(r) is understood as the bi- 
variate operator e(r|r') = e(r)8^ d \r — r') so that (eE)(r) = 
/ d<V e(r|r')E(r'), Eq. <[3J) reads: 



E = E b + G l/d {8e/e b )E. 



(3.3) 



From this point on, we depart from the usual treatment. A 
modified Green operator Gg is considered, in which the sin- 
gularity at the origin is "renormalized" by introducing [|3"Tll an 
arbitrary "inner" depolarization parameter £, of a more funda- 
mental nature than £ e . By definition, 

G 1/d (r) = -*l<5(r) + G,(r). (3.4) 

Writing the prefactor of the local (Dirac) part of Gg as 

A = £-l/d, (3.5) 

one has by Eq. d3.21 i and definition (13.41 , 

Gg(r) = A\S(r) - lim 6{r - ^ — -rO - drr). (3.6) 

>;->o b d r a 



E/ 



and a modified permittivity contrast 



u e (r) 



Mr)/e b ] 



E, 



e(r) - e b 



(3.7) 



1 + e[6e(r)/e b ] fe(r) + (1 - t)e b 
Eq. ( 13.1b in transformed into the equivalent equation ll3lll 



(3.8) 



El — En 



Gg ug Eg 



(3.9) 



For I = 0, uo(r) is the scaled dielectric contrast. For £ = 
Wi/d is (up to a prefactor vEb) the dielectric polarizability 
of a spherical cell, and Ei/d(r) is the local (Lorentz) field 
impinging on it 0]. For £ 7^ 1/d, some screening effects 
are implemented at the level of the cell. This possibilit y of 
splitting the Green function has been noticed previously [47] 
to the purpose of improving convergence in solving Eq. (13.1) 
by iterations (see also [4] p. 302). Although related, our aim 
is somewhat different. 



B. Multiple-scattering expansions 

Both ( 13.3b and ( 13.91 ) are continuum analogues of scattering 
equations for finite-size scatterers. It proves convenient to em- 
phasize the connection by casting the continuum problem into 
a multiple-scattering framework. The scattering potential ug y 
of the material element at y is introduced as 

u, y (r|r') = v Ui (y)6(r - y)S(r' - y) I. (3.10) 

Setting Gg(r\r') = Gg(r — r') Equ. ( 13.91 ) is rewritten as 

Eg(r) = E (r) + / -^(fV 2 (fVi G^(r|r 2 )u £y (r 2 |r 1 )E^(r 1 ). 

(3.11) 

In order to distinguish between integrations over 'in' and ' out' 
space variables ri 2 and summations over scattering cells, let 
ugy = ugj and make the formal replacement J(d d y/v) — > 
J2j- Then (O takes the form 



Eg — Eq + Gg ugj Eg. 
3 



(3.12) 



Define now the Green function G associated to E( by the 
equation Eg = GGJ 1 Eq, where G~[ x En is the source. Then 

G = Gg + Gg^ugjG. (3.13) 



The individual transition operator, also called T-matrix (see 
l42ll for a comprehensive list of references), completely char- 
acterizes the polarizability properties of the i th scatterer. It 
reads 



U = - Ggugi) x . 



(3.14) 
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Introducing the Green function for the local field impinging 
on the scatterer, namely, G; = (1 — GfU^G; the polar- 
ization of the scatterer Pi = tiG L ; and finally its "dressed" 
T-matrix, Tj = PiGj 1 , which accounts for corrections due 
to the other scatterers, one obtains from (13.131 the familiar 
multiple-scattering equations i37ll : 



G — Gi + GffGi 1 



(3.15a) 
(3.15b) 



where T is the T-matrix of the whole system. We observe that 
since scattering events occur in succession between different 
scatterers, the Dirac term of Gg is irrelevant in T. Thus, we 
replace Gg by Gi/d in (13. 15bb to indicate that this propagator 
does not depend on I. 

Next, a n-body ("virial") expansion of the form T = 
Sn>i T*-™) is written down, where the n-body term 



fin) = 



1 / rp(n) 



(3.16) 



is a sum over all possible scattering sequences on n dis- 
tinct scatterers. The prime indicates that the sum is carried 
out first over label i-y, then over %2 ^ i\, etc., then over 
in & {h, «2, • • • , in-i}, so that depends on la- 

bel ordering. This question was solved long ago by Peterson 
and Strom [48]. We provide in Appendix [A] another different 
approach to their result. Specifically, we show that: 



T. 



(n) 



(ii,i2,. ■■,in ) ^(ii^T-'-in-i)^^-! 'd^iri 

x (l — Gi/dSi 1 ...i n _ 1 Gi/dti n ) 
x + ,. in _ 1 ), (3.17a) 

Si!. ..i n = Gy d [(l + Gl/ d t in ) 

x (l — Gi/dSi 1 ...i n _ 1 Gi/dti n ) 
x{l + G 1/d S il ... in _ 1 )-l], (3.17b) 



(!) - O. - 



(3.17c) 



In particular, the two-body term, to be used hereafter, is 

IH: 

= thGi/dtiz 

+tiiGi/dti 2 Gi/dti 1 Gi/dti 2 (l — Gi/dt^Gi/dtiz) 

+ thGi/dUzGi/dth (l — Gl/dtisGl/dtii) ■ (3.18) 

In the last writing, the term G\/dti 2 has been singled out for 
convenience, in the perspective of using Eq. (I3.25bb below. 
Although we disregard it for simplicity, the three-body term 
ll48ll has been considered by Cichocki and Felderhof lIBoll . 



C. Ensemble averages, self -energy and effective permittivity 

To address substitutional or positional randomness, ensem- 
ble averages over disorder, denoted by (■), are carried out 



113 811 . Due to statistical homogeneity, averaged operators are 
translation-invariant. Introducing the complete scattering po- 
tential Ug = uu associated to the whole set of hetero- 
geneities, the so-called self-energy (or coherent potential) op- 
erator Yig of the averaged Green function (G) Il44ll is defined 
by (U(G) = Y,g (G), and from (3. 13) follows the Dyson equa- 
tion jH 



(G) = (GJ 1 - 



(3.19) 



The kernel E^rjr') = E^(r — r') possesses one local and one 
non-local part. It can be written as the Green function (13.2b in 
the form 



E,(r) = Ef°\8(r) 



v-nloc 
1-t 



(3.20) 



where E)? c is a scalar, and where the same principal-value 
prescription as in (13.21 i applies to the non-local term. The ne- 
cessity of distinguishing between the local and non-local part 
will show up when comparing theory to our reference data. 

Similarly, the effective permittivity e e is a non-local opera- 
tor 115 211 . of the same generic form, and is defined through the 
equality (eE) = e e (E); that is, 

(e(r)E(r)) = J dV e e (r - r')(E(r')). (3.21) 



Definitions of Eg, e e and T,g lead to 113 111 

e e = e b [l + (1 - *)£/][l - ffi/]- 1 . (3.22) 

This formal equation becomes algebraic when Fourier trans- 
forms of the kernels are used. In particular, it holds for the 
volume integrals of the kernels. We recall that the effective 
£ e (/) is obtained from e e by means of ( 12.91 ). 
The configurational average of ( 13. 15ab yields 

<G) = G e + Gg(T)Gg. (3.23) 

Comparing this equation with (13 . 1 91 > provides the relationship: 



H e = (T)[l + G e (T)Y 



(3.24) 



For an infinite system, this equation is merely formal and only 
has a meaning as a perturbative series, because (T) involves 
conditionally-convergent integrals. Expanding (13.24b . and us- 
ing ( 13.161 ) and d3.17| ). gives in the form of a n-body expan- 
sion, where one- and two-body contributions are read from 
(T3~T7cb and (EH: 



(n) 



where E 



(i) _ 



(3.25a) 



w= E( T (5 2 ))-E<*«i) G ^^); etc - < 3 - 25b ) 



(2) 

The long-range part of the second term of E) ' cancels out the 
first term of the average of (13 . 1 8b . The remaining terms in- 
volve only absolutely-convergent integrals, with an integrand 



decaying at least as G x 



Id 



r 2d as r 



This prop- 



erty holds at each order (the perturbative expansion of E has 
a special - ordered - type of cumulant structure ll53l UU), 
which implies that E is independent of the sample shape in 
the infinite-volume limit 15415511 . 
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D. Parameters s and q, and clustering 

Before computing the self-energy, we remove some inde- 
terminacies of the continuum theory. The latter admits two 
natural adjustable parameters. A first parameter s > stems 
from remarking that 77 in Eq. (13.21 i must be of order a, which 
we write r\ = sa. We interpret s/2 as a measure of an effective 
radius of neighboring inclusions, to be further constrained in 
Sec. IIV Al The value s = 2 (77 = 2a) corresponds to set- 
ting a hard-sphere-type exclusion distance between polariz- 
able point inclusions. However, we allow here for smaller or 
larger values to tune the strength of dipolar interactions, so as 
to compensate for the lack of explicit higher-order multipoles 
in interactions between finite-size inclusions 156, ,5711. 

To work out the above equations, we need an expression 
of ti, defined formally by Eq. d3.14t . Dropping the index i 
for brevity, and expanding, one has t = u^'^2 lk>Q {GiUg) k . 
Powers are evaluated by means of definitions (13.6b and (13. 10I ). 
but this involves squares of Dirac distributions: consider for 
instance ui{GiUf) 2 , which reads 

f 4 

/ ]^d^ fe U£ y (r|ri)G^(ri - r 2 )u£ y (r2|r 3 ) x 
3 k=i 

x G e (r 3 - r 4 )u£ y (r4|r') 
= v 3 u e (yfGi(0) 2 S(r-y)S(r> -y) 
= v u e (y) 3 A 2 [vS(0)] 2 6(v - y)8(r' - y). (3.26) 

As discussed in Appendix|B] we use the prescription 

vS{r) 2 = qS(r). (3.27) 

The number q > 0, a mathematical and physical necessary 
addition when A 7^ 0, is the second parameter of the theory. 
We can then write vS(0) = q. The "infinitely large" number 
v~ x gives the physical "order of magnitude of (5(0)" in this 
problem. Thus, 

t = u e J2(l A ) ku e = M 1 - q^u,,)- 1 = uj, (3.28) 
where 

£ = £-qA = (l-q)£ + q(l/d); (3.29) 
that is, with t(y) = uj(y) 

t y (r|r') =vt(y)6(v - y)S(v' - y) I. (3.30) 

We assume that < q < 1, so that I can be interpreted as a 
weighted average of £ and 1/d. 

We propose the following interpretation of g, illustrated by 
Fig. 12 where ellipsoidal shapes are meant to indicate that in- 
clusions have an effective depolarization factor different from 
1 jd (the orientation of the ellipsoids in the drawing is irrele- 
vant). Each scatterer ti is viewed as a aggregate of screened 
polarizable elements of polarizability proportional to ui, of 
same permittivity, embedded in medium Eb- Their degree of 
clustering is adjusted through q. We interpret the latter as a 
coverage/spreading parameter for aggregated elements. In the 




•<"*! Low coverage: •••"••C-v High coverage: 

FIG. 2. (Color online) Clustering interpretation of the s and q pa- 
rameters. 



figure, the small bar associated to A represents some spread- 
ing of <5(r) (see Appendix IB1. When q = 0, t = ug and the 
elements are considered separately; when q = 1 they gather 
as a compact isotropic (spherical) inclusion of polarizability 
t = ui. Setting £ = 1/d suppresses the influence of q: 
the aggregate reduces to one single spherical inclusion in this 
case also. The aggregate is represented as a whole in a mean- 
field way in the sense that electrostatic interactions between 
its components are approximated by the local part A of Gi . 
Moreover, aggregates are treated as point-like polarizable ob- 
jects [Eq. ( I3.30H when it comes to considering their mutual 
dipolar interactions via Gi or Gi/d (of which only the non- 
local part is relevant here; see Sec. lIHBl . 

To summarize, introducing < q < 1 makes the above 
theory a simplified one of interacting aggregates, in which the 
difference between £ and £ distinguishes between inclusions 
and aggregates thereof. 



E. Self -energy to two-body order 

The self-energy follows from (13.251 . Considering only 
volume-integrated kernels (with a slight abuse of notation), 

E (i) _ E (i)ioc = j d<V£W(r-r')=y d<V^](^(r|r')) 
= j d<VdV(r - yW - y)(t) - (t). (3.31) 

Likewise, from GTiM and d3~25bl , E< 2 ) = £( 2 ) loc + £( 2 ) nloc 
with 

s(2 )loc = _ A(<) 2 + U I ^ G l/» \ (332a) 

/ v 2 t 2 t 2 G : * (r) \ 

E (2)„loc = U I 1/dV ) \ 

Jr> v \l-vHtG 2 1/d (r)/' 

where rj — sa, and where r is the separation vector between 
two statistically uncorrelated cells of volume v. We denote 
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their polarizabilities by t and t to distinguish them. The con- 
figurational average over permittivities must be carried out in- 
dependently on these quantities. 

The first two terms of (I3.18l l are built on scattering se- 
quences that start and end on the different scatterers i\ ^ i 2 , 
and thus contribute to E( 2 ) nloc , which stands as a non-local 
susceptibility. Instead, the third one is made of closed scatter- 
ing sequences that start and end on the same scatterer i\, and 
so contributes to the local part I]( 2 ) loc as a renormalization 
("dressing") of the T-matrix of individual scatterers 114911 . In- 
deed, for A = 0, in the diagrammatic representation ll2ll l40ll 
where a line stands for a "propagator" G\u and a 2n-legged 
dot stands for the nth power of a T-matrix t, 



s (l)loc = A 



£j(2)loc 
^(2)nloc 




(3.33a) 

(3.33b) 
(3.33c) 



As recalled in Sec. lHICl the absolute convergence of the inte- 
grals in £( 2 ) results from the fact that the latter involves only 
terms with at least two propagators. When A ^ 0, the addi- 
tional contribution — A(i) 2 , which stems from the local part 
of Ge in the last term of (13.25bt . is counted within E (2)loc 
rather than £( 2 ) nloc , in view of definition ( 13.201 i of the local 
and nonlocal parts of the E operator. 

In spite of the shorthand scalar notation employed in (13.321 
(all terms commute) matrix inverses are required. They follow 
from writing the nonlocal part of G\/d as a linear combination 
of the projectors f f and I — f f . Carrying out the integrals over 
r in d332|> yields ifssh 



E (2)ioc = _ A ^2 + d -2( d _ / t (tt)?g(z)) ,(3.34a) 

s (2)nloc = d -2( d _ y (tt h ( z ^ ) (3.34b) 

g(z) = tanh. _1 (z) + tanh _1 ((d- l)z), (3.34c) 

h(z) = ~l 0i 



1-z 2 



2 b 1- (d- l) 2 z 2 ' 
z = Vtt/(ds d ). 

We note for further use that, with principal determinations 

l + z 



(3.34d) 
(3.34e) 



g(z) + h(z) = log 



1 - (d- l)z' 



(3.35) 



Due to scale invariance, size a is absent from these equations. 

Gathering one- and two-body terms, the explicit expression 
of E^ after having taken averages is as follows. We set p\ — 
(1 — /), P2 = f, and introduce t\ and such that 



£1,2 - £b 



1,2 



&1.2 + (1 - tyb 



(3.36) 



The quantities t and t in Eq. (13.321 > take on values t\ or £2 
with respective probabilities pi and p2- Introduce moreover 

z n = t 1 /{ds d ), z 12 = y/~ht2~/{ds d ), z 22 = t 2 /{ds d ). Then, 



from Eqs. (ED) and (l3~34l with (t) = pih + p 2 t 2 , 

£j? c = (t) - A(t) 2 + d- 2 (d - 1) [(p^fgizn) 

+ {P2t2f g{z22) +PlP 2 (*l +t 2 )(*l*2)'5(2l 2 )], (3.37a) 

E? Ioc = d- a (d-l)[(pit 1 ) a /i(z 1 i) 



+ (p2t 2 ) 2 Hz 22 ) + 2( Pl t 1 )(p 2 t 2 )h{z 12 )] . 



(3.37b) 



We remark that h(z) = for d — 2, so that in two dimensions 

— Vloc 

£ = 2j £ . 

Two-body interactions terms in the effective permittivity of 
composites have been considered by numerous authors, many 
of who focused on interactions between spherical inclusions 
of finite size (HIH. Here ' takin 8 d = 3, £ = 1/3, s = 2, 
£(, = E\, and ti — a/(e\v) where a is a polarizability, and 
t 2 = 0, and expanding e e to order 0(f 2 ) returns a known ex- 
pression of the two-body correction in the effective permittiv- 
ity of a suspension of polarizable point inclusions distributed 
according to the law of a hard-sphere gas ll57[l60ll . An often- 
cited expression for this term Irolll is only an approximate one. 



IV. EFFECTIVE-MEDIUM CONDITIONS 

In this section, the general theory is completed by s.c. con- 
ditions, and exploited. Several schemes are possible. For clar- 
ity, a non-self-consistent setting is considered first. 



A. Theory of the Clausius-Mossoti type 

When used with the volume integrals of the kernels, Eq . 
(13.221 1 resembles the Clausius-Mossoti (CM) formula JH], 
in which the polarizability of inclusions replaces the self- 
energy, and constitutes a generalization of the Maxwell- 
Garnett effective-medium formula [63]. The latter holds for 
dilute systems of spherical inclusions embedded in a matrix. 
It is retrieved by letting 1 = 1 /d, fixing the background per- 
mittivity to that of the matrix, and keeping only in E the 
one-body contribution (13.25ab . Similar many-body general- 
izations of the CM formula by means of cluster expansions 
have previously been worked out by Felderhof and co-workers 
among others. 

Whenever E^ 7^ in Eq. ( 13.221 ). which produces CM-type 
estimates, a distinction must be made between the backgound 
"inner" permittivity Eb and the overall one e e . It follows that 
the "inner" £ of the theory needs to be distinguished from the 
"effective" one, £ e (f) obtained from e e by means of Eq. ( 12.91 ). 

In the rest of this Section and in the next one, we fix £ to its 
usual value 1/d. This eliminates altogether the q parameter, 
see Sec. IIII Dl and implies that A = and 



d(si,2 - e&) 
£1,2 + {d- !)£&' 



(4.1) 



Dilute behavior in the CM-type approach is then as follows. 
Setting Sb = £1 the 0(f) term of £ e is readily obtained. Using 
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identity ( 13.35b . the result reads 



e? M (f) = - d ' /_] 



<2 



d 2 
£2 ~ £l 



d( 2 - log 



1 + C 2 /s a 



l-(rf-l)C 2 /s d 



£2 + (d - l)e% ' 



(4.2) 



The medium being symmetric, the 0(f — 1) behavior stems 
from interchanging e\ and £2, and from replacing / by 1 — /. 
Letting 4(/) - l/d + a f + O(f 2 ) = 1/d + <n (/ - 1) + 
0((f — l) 2 ), the dilute-limit slopes of the £ e (f) graph in the 
infinite-contrast limit £2 £1 follow as 



.CM 



.CM 



d-l 

d-l 



d - lo S 3T 



1 



s" - (d - 1) 
d , s 



d-l 



loe 



d 1 1 



s d - (d-l)" 1 



(4.3a) 
(4.3b) 



By self-duality U, a% M = of M for d = 2. 

The cut of the logarithm in ( 14.21 i materializes the two-body 
resonance spectrum [Hill of the theory, and such resonances 
should not occur for positive e\ t i\ that is, £ e must be real in 
statics. Assuming that s does not depend on volume fractions, 
the requirement that the infinite-contrast limits of cto.i be real 
implies the constraint 



s > s n 



(d-l) 



l/d 



(4.4) 



This is a direct consequence of requiring that (d — 
1) max[|zn|, \z\2\, \z22W < 1 in Eqs. (13.371 1 to avoid the cuts 
whatever d and the contrast. 

Fig. |3la) represents the functions a™( s )- Slope 00 blows- 
up logarithmically near s m ; n . Thus, for d = 2, negative slopes 
cto.i such as in Figs. [Tfa) and (c) are possible for s close 
enough to 1. For d = 3, negative slopes ctq are obtained for s 
close to s m ; n (d = 3) ~ 1.25992, and the theory admits pos- 
itive (7 1 slopes: this is evocative of the PDA behavior of Fig. 
|TJd) in the dilute limit, but the corresponding <j\ values, of or- 
der 0.2, outstrip those of simulations. Besides, this approach 
cannot reproduce the weak negative slope at / = 1 of Fig. 
Eb). 




o- 1 (d=3) 
CM-type theory 



1.0 1.5 2.0 2.5 3.0 3.5 4.0 

s 




FIG. 3. (Color online) Slopes 00,1 in (a) the CM-type theory, and (b) 
generalized self-consistent theory, with £ — l/d, vs. parameter s. 
Vertical dashed lines represent the s m m lower bounds. 



B. Generalized self -consistency with £ — l/d. 
Parameters a and /3. 

We keep using £ = l/d. The exceeding high slopes ( 14.31 ) 
are a consequence of their going to a finite limit as s — s- 00, 
because of the term d^2 within the brackets of Eq. ( 14.2b . This 
term is traced to the one-body contribution X^ 1 ) to the self- 
energy, i.e., it already exists in the standard Maxwell-Garnett 
formula interpreted in terms of £ e . Powers of £W are gener- 
ated to all orders of perturbations when expanding Eq. (13.221) . 
The self-energy expansion d3.25t is essentially an asymptotic 
one, and the CM-like expression (13.221 ) does not perform well 
in reorganizing it into a physically-meaningful effective per- 
mittivity for all volume fractions. 

A widely-used reshuffling device is to determine £& self- 
consistently by the CPA condition £ = + E (2) = 0, 
which implies that e e = £&. However this is not the only pos- 
sibility: we can also consider the restricted one-body version 
£W = ll43ll . i.e., the BL condition; or even an intermediate 
"local" condition S loc = + £( 2 ) loc = 0, etc. To handle 
them all, and much more, we introduce new parameters a and 
j3, and put forward the generalized s.c. condition 

+ aE (2)loc + ^£( 2 )" loc = < a, /3 < 1. (4.5) 

Parameter f3 is irrelevant for d = 2, according to our remark 
following Eq. ( I3.37bb . With the above condition, and setting 
ST, = (1 - a)£( 2 ) loc + (1 - /3)£( 2 > nloc , Eq. (|3~22| | reduces to 



£b 



1 + (1 - l/d)ST 
l-ST/d ' 



(4.6) 



which, unless a — (3 — 1, remains of the CM type in spite of 
the s.c. condition. 

As far as the O(f) term in £ e is concerned, the outcome 
of condition ( 14.51 ) is independent of (a, (3); this property does 
not hold for £. With (2 as m ( 14.21 ), working out the dilute 
expansion indeed gives 



v» - 1 



1 d-l 



d 2 



4 



log 



1 + C2/S 



i-(d-i)c 2 / S fl 



(4.7) 

It is important to remark that with self-consistency, the slope 
at / = becomes a function of (2 / s d , which turns parameter 
s into a polarizability rescaling factor. Infinite-contrast slopes 
follow as: 



d-l 

d-l 

d 2 



d 1 



s a + 



1 



(d- 



ds- d 

Jd~T) 



- log- 



(d-1)" 



(4.8a) 



(4.8b) 



They are drawn in Fig.|3b). They vanish in the limit s^oo 
where two-body interactions are suppressed, so that the slopes 
are now a pure two-body effect. However, whereas o\ takes 
on more realistic values, it still cannot be made negative for 
d = 3. 

Percolation occurs through e^. A second-degree polyno- 
mial equation for the percolation threshold f c is obtained as 
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FIG. 4. (Color online) Dimension d = 3. Dependence on ft (as indicated) of £ e (f) in the generalized s.c. approach at infinite contrast, with 
s = s m in + 0.05 and a = 0.00, 0.33, 0.67 and 1.00 in each figure (top to bottom). 



the infinite-contrast limit of (14.51 1. letting Ei = first, then 
Eb = 0. For brevity, we no not reproduce its lengthy coeffi- 
cients. The solution for d = 2is/ c = l/2 irrespective of the 
parameters. For d — 3, the model percolates at the BL thresh- 
old f c = 1/3 if a = f3 = 0; otherwise, the leading-order term 
of the equation for s near to s m i n provides the asymptotic es- 
timate 



fc 



6 + alog5-/31og(5/4) 



-l 1/2 



(a + p)\log(s- 



(4.9) 



The lowest value of f c is obtained with a = (3 = 1, at fixed s, 
and f c can be made as small as needed by letting s —> s min - 

Fig. E] illustrates typical consequences on £ e (f) of the gen- 
eralized s.c. condition, in the three-dimensional case where 
variations are most conspicuous. As long as self-consistency 
involves two-body interactions, the threshold is lowered with 
respect to the BL value, with an intricate dependence on a and 
f3. Their influence on l e is confined to a definite region around 
f c , because the dilute-limit slopes do not depend on a and (3. 
The effect of removing part of £( 2 ) nloc from the s.c. condition 
is interesting: in Fig. Ufa), a graph shape with an upward cusp 
at f c , akin to that of Fig. |TJd) (highest contrast), is produced 
with (3 = and a small amount of £( 2 ) loc ( Q = 0.33); on 
the other hand, a downward cusp results from injecting a high 
amount of £( 2 ) nloc [fj = 1, Fig.gtd)]. 



C. Generalized self-consistency with variable I 

We can now discuss the effect of the "inner" depolarization 
variable £. Relaxing the condition £ = 1/d, the s.c. condition 
(14. 5t for Eb still holds, but this time expressed in terms of 
with £ arbitrary. Since A = I — 1/d 7^ 0, the self-energy 
has one additional term in its "loc" part (13.37a| i, and the q 
parameter introduced in Sec. IIIID I becomes operative. The 
theory of Sec. lIVBl is retrieved if q = 1. To determine £ we 
impose the supplementary s.c. condition, 



0, 



(4.10) 



which, in the interpretation of Sec. MI Dl makes E\> the local 
effective medium surrounding the elementary screened ele- 
ments within aggregates. This is Eq. (II . lb . with e e and £ e 



replaced by Eb and £, respectively. Its solutions are 



1-1 
1 



■ ± y/e{i - £)T X + 



T=^[{l-f-£) + [f- £)xl X = ei/ei. (4.11) 

Solutions £ follow from using these expressions in the gener- 
alized CPA condition d4.5t . If a = 1 (case d = 2 where j3 is 
irrelevant), or a — (3 = 1 (case d > 3), then E = so that 

£ e = Eb and £ e = £. 




(b) s 




physical branch 






s=3 




q=1.02 10~ 3 




o=1.0 



0.0 0.2 0.4 0.6 0.8 1.0 
f 



0.0 0.2 0.4 0.6 0.8 1.0 
f 



FIG. 5. Dimension d = 2 and contrast 10 s . Typical conformation of 
the branches of I for d = 2 with (a) a 'bad' case; (b) a 'good' one 
(see text). 

Depending on parameters {/, s, q, a, /?} the theory admits 
up to three real solutions. Admissible values of {s, q, a, (3} 
must be such that a continuous real function £(f) exists for 
all contrasts in the interval / £ (0,1), with endpoint values 
i = 1/d. We call it the "physical branch". A criterion that 
generalizes (14. 41 > to the whole parametric domain seems out 
of reach. Still, by requiring real solutions for positive permit- 
tivities and any < I < 1, an argument similar to the one 
used in Sec. lIV Al leads to the constraint 

qs d >(d-l), (4.12) 

which extends d4.4| > to values of q ^ 1. The actual permitted 
domain depends on (a, (3) and is certainly wider, since solu- 
tions I take on values in a much more restricted interval. Thus, 
the above constraint is too restrictive in practice. 

Numerical experimentations show that for d — 3, percola- 
tion thresholds markedly lower than 1/3 are obtained for pa- 
rameter values close to the boundary of the domain, much as 
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in the simpler case of the previous Section, illustrated by ex- 
pression ( 14. 9i of the threshold. As a rule, close to boundary, 
the smaller q, the larger s, which is consistent with Eq. (14.121 
and with the interpretation of Fig. [2] the more dispersed the 
clusters, the larger their effective radius. But for lack of any 
simple analytical expression of the boundary of the allowed 
domain, no systematic study of the threshold in d = 3 was 
made. 

Fig. [5] displays some typical conformation patterns of so- 
lution branches of £ for d = 2. Parameters, as indicated (re- 
call that f3 is irrelevant for d = 2), were chosen such that 
£'(f = 1/2) ~ 1. The expression of £' (/ = 1/2) is easily 
deduced by series expansions. The solid (resp. dashed) lines 
represent solutions arising from using the 'plus' (resp. 'mi- 
nus') branch in Eq. ( 14.1 It . Lines interrupt themselves where 
solutions become imaginary. Bad parameter values, such as 
in Fig. |5ja), lead to discontinuous real solutions £(f). Fig. 
|3b) illustrates a "good' case, with the physical branch indi- 
cated. Although lying close to / in a region around f c = 1/2, 
the function £(f) could not be made close enough to repro- 
duce the double-threshold effect of checkerboards mentioned 
in Sec. Ill Al The physical branch was always found to be gen- 
erated from the 'plus' solution of (14.1 It in all cases exam- 
ined where it exists. Whereas in the infinite-contrast limit and 
for some special parameter sets, joining 'plus' and 'minus' 
branches can produce a composite continuous 1(f) graph with 
endpoints 1/d, such special solutions do not survive at lower 
contrasts, and cannot be considered physical in the present 
context (although the situation might change upon including 
three-body interactions). 

With (,2 as in (14. 2b , and introducing 

a g = [q + (l-q)a}/q, P g = [q + (1 - q)/3]/q, (4.13) 
and the function (z) = dz—g(z), the dilute expansion reads 



C'(/) 



1 d- 1 



d 2 



a i92 [ — d 



/■ (4.14) 



When q = 1, it reduces to ( 14.7b thanks to identity (13.35b . The 
infinite-contrast slopes at / = 0, 1 become: 



d 2 

d-1 
~d?~ 



d-1 



d-1 



(4.15a) 



(4.15b) 



With q 1, the slopes now depend on a and (3, with an over- 
all scaling by 1/q. The effect of varying a q and (3 q can be 
studied by factoring out j q = (a 2 + (3 2 ) 1 / 2 , and by introduc- 
ing the angle < <fi < n/2 such that <fi — arctan(/3 9 /a 9 ), for 
various values of t/>. 

In two dimensions the common slope is a q times that of 
Fig-Ob). Fig.[6]displays graphs of the scaled slopes <JQ C f/^f q 
vs. s in three dimensions for some values of </>. The important 
point is that values of u\ with either sign are now available. 




0.0 0.2 0.4 0.6 0.8 1.C 

S — Smin 



FIG. 6. (Color online) Dimension d — 3. Scaled slopes vs. parameter 
s, in s.c. theory with variable £, for 20/tt = 0.0, 1/3, 1/2, 2/3, and 
1. Dashed: cto; solid: ai. Arrows indicate the direction of variation 
as <j) increases (color online). 



When s d 3> 1 and q <C 1 the slopes behave asymptotically as 



jf*~ -a/{6qs 6 ) if d = 2, (4.16a) 

-P/(3qs 6 ) if d = 3,/3^0 

-2a/(3qs 9 ) if d = 3,/3 = 0, { ' 

l3/(12qs 6 ) if d = 3,/3^0 

-a/{12qs 9 ) if d = 3,13 = 0. { ' 



Figure [7] compares the simulation data of Fig. Q]with theo- 
retical plots drawn using suitable parameter values in the s.c. 
equations. Parameter values were constrained by imposing the 
theoretical slopes at / = 0, 1 to match on average those read 
from the data, and by requiring a high slope at / = 1 /2 in the 
two-dimensional case, which is achieved by taking s arbitrary, 
but large. This does not uniquely determine the parameters, so 
some variations are admissible. 

In both two-dimensional cases of Figs. [Tfa) and (c), good 
agreement with the data is obtained for low to moderate con- 
trasts, but the high-contrast data values for / close to f c can- 
not be matched. Actually, the peak values of £ e (f) at high 
contrast could have been reached by using other parameters, 
but at wrong volume fractions and at the expense of the low- 
contrast fits. It can be shown that for d = 2 the theoretical 
slope at / = 1/2 is at most a(f = 1/2) = -4<r(/ = 0, 1), 
which is obtained for s large. This explains why the high- 
contrast data cannot be matched. Indeed, a(f — 1/2) should 
ideally be equal to 1, whereas a(f = 0, 1) ~ —0.05 in both 
sets of two-dimensional data. 

The best overall fit is obtained in the d = 3 RRN case, 
Fig. 0b). Good match is obtained in the three-dimensional 
PDA case of Fig.[7|d) also, but in the low -permittivity region 
(/ < fc) only. In the high-permittivity region (/ > f c ), the 
dilute slope is well reproduced, and a cusp near /+ at high 
contrast is retrieved; however, £ e {f) is too low, and cannot 
be turned into matching the data by varying parameters. This 
case requires a q value much larger than for the other three 
and, correspondingly, a relatively small s value. 
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FIG. 7. (Color online) Theoretical plots of £ e (f) obtained from the s.c. theory with variable I (solid) fitted on the data of Fig. Q] (dots). 
Contrasts and parameters values are as indicated. 



V. CONCLUDING DISCUSSION 

We summarize our findings, and discuss possible improve- 
ments and directions for further work. First, we proposed 
to interpret effective-permittivity data in terms of an overall 
depolarization-coefficient function £ e , which proved a sensi- 
tive probe to highlight differences between models in a neat 
way, for all concentrations and contrasts. 

Our theory rests on the introduction of an inner free depo- 
larization coefficient I, Together with the usual free back- 
ground permittivity, e^, these quantities are determined by 
coupled s.c. equations: a Bruggeman-like equation on the po- 
larizatibility of clusters, and a CPA-like equation on the over- 
all self-energy. They reduces to Bruggeman's with I = 1/d in 
the absence of many-body corrections. The possibility of in- 
troducing £ and stems from an overall invariance property 
of multiple-scattering theory under such parametrizations. Al- 
though we assumed et, and £ to be scalars, it is clear that the 
most general reparametrization should involve tensors. This 
possibility was not considered, but would be necessary, e.g., 
to compute field fluctuations, since this involves consider- 
ing anisotropic perturbations to the permittivities [66 j . As it 
stands, the theory only applies to cases where the one-body 
term in the effective permittivity is isotropic in the dilute limit, 
given by Eq. ( 12.2b . 

Two empirical parameters s and q were introduced to de- 
scribe clusters. While values £ ^ 1/d suggests that the rel- 
evant clusters of the theory are not spherical, computations 
with two-body terms were carried out with a spherical exclu- 
sion volume (parametrized by s), for simplicity; we note that 
a parameter n, similar to our s, was previously introduced by 
Cichocki and Felderhof in two-body integrals to the purpose 
of studying scaling relationships in corrections to the CM for- 
mula J6l. 

Allowing for £ ^ 1/d involves products of Dirac func- 
tions. The latter arise from a rough treatment of clusters as 
point-like polarizable entities, within which dipolar interac- 
tions occur only through the local part of the Green function. 
The resulting mathematical ambiguity was handled by intro- 
ducing q, interpreted as a covering parameter between inner 
inclusions. The main operational role of parameters s and q is 
to modify the three-dimensional percolation threshold of the 
theory. Inasmuch as they represent microstructutral features 



of the medium, these parameters are akin to the geometrical 
ones introduced by Miller [25], although we cannot claim any 
precise connection at this point. 

We moreover needed to introduce two supplementary pa- 
rameters a and (3 to handle independently the local and non- 
local two-body parts in the CPA self-consistency condition on 
the self-energy. When only part of the latter is canceled out, 
the theory stands as intermediate between CM-like and BL- 
like approaches. This additional flexibility was required to 
produce realistic £ e (f) functions. 

Although our emphasis was on versatility, the present the- 
ory could be rigorously reformulated in the discrete frame- 
work adapted to RRNs, for which the two-body interaction 
term is known exactly (without needing to introduce param- 
eters s and q) 111 U l49l 15811 . in order to investigate precisely 
the role of parameters a and (3. We also remark that while 
the percolation threshold can be adjusted for d = 3, where 
it depends on all four parameters, it remains stuck to its ex- 
act bond value f c = 1/2 in d = 2. This does dot allow one 
to handle site percolation [5]. Improving the behavior near 
f c , as well as attempting to reproduce site-percolation effects 
within the present theory, would presumably require including 
three -body interactions Il50ll . It would be possible, if needed, 
to introduce additional parameters similar to a and f3 when 
including their contribution in the generalized CPA condition. 

Finally, the multiple-scattering formulation in the contin- 
uum allows in principle for extensions to elasticity Il68ll . or to 
the frequency domain, including wave propagation dynamics 



Appendix A: iV-body expansion in multiple-scattering 
framework 

To prove (13.171 i. we first demonstrate that 

n _j 

T (h,i2,-,in) = T (h^l-,in-l) GtUU - i 1 " ^ U «p) ' 

p=l 

(Al) 

for n > 2, the recursion being initiated with T$ = U. We 

start from the expression T = J2i u n{^ ~ @t Ylj 1 > or> 
tained from the equivalence between (I3.13l l and ( 13. 15ab . and 
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to which we apply systematic transformations. Introduce 

di = (1 - Gzuu)~ l , Aj= (l — Gf, uikOLj) . (A2) 



Then, with =U = uucn 



T = uu [ 1 — Giuu — Gi 22 uij i j = ^ i* ^4 



j#i 



EM 1 + Gi uejdiAj 



E T S )+ E4) )G ^ a ^- 



(A3) 



Let now ay = [1 — Gt(ua + . One has 

A^ 1 = 1- Gputj a,i- Gi 2_j u ^ a i 



l-Gt^ u lk a. L (1 - GtUijCti) 1 (1 - Gluten) 



i - Gi y^ ".</.",, )",,'",. 



(A4) 



Whence, with A y = (l - G f Z)*^ w « a «j ) > and 



r (?j) ^T^GiUejUij, 



(A5) 



it follows that 

= E T W + E T S) E 7(8)^ W A« , 

i * i 

(A6) 

The progression from (1A3) to ( |A6) represents one transfor- 
mation step. Going on by applying to the last term of JA61 > a 
transformation similar to (1A4) . namely, 

Ar} = 1 - G^Ufc dij - Gi 22 u " a n '» ( A7 ) 
and so on, proves dAl) . Proving ( 13.17) is now simple. Assume 



(13. 17ab to hold at rank n, and evaluate Tt 



<n+l) 
(ilia— in+i) 



. Let 



n _, 



P =i 



Because of dAll ), it suffices to show that 

Bn+i (A8) 
= ^n+i (l — GiSi 1 i 2 ,,,i n Giti n+1 ^j (1 + GeSi 1 i 2 ,,,i n ) 



with Si 1 i 2 ...i n expressed in terms of &£ 1 .£ 2 ...j„_ 1 as in d3.17b) . 
knowing that 

Uii„B n = ti n (l — GeSi 1 i 2 ...i n _ 1 Giti n ) 

x {l + G t S lll2 . ..*„_,). (A9) 



Since by definition [cf. ( I3.14H u« n = £i„(l + G^,- n ) 1 , 
implies 



-B n 1 — (l + G£S' ili2 ... lii _ 1 ) 

x (1 - G f S lll2 ... ln _ 1 G e U n ) (1 + G^J- 1 . (A10) 
Therefore 

B n+i = B n X - GeU n+1 (l + Gi>t ln+i y } 
= (1 + GiSi^i 2 ...i n ) 



x [1 - (1 + G e S lll2 ... ln )G e t ln+1 (l + Gzt tn+I y 1 } 

= (1 + G^S'i 1 i 2 ...i ii ) 
x (l-GiS ili2 ... in G e t in+1 ) {l + G e t in+1 )-\ (All) 



where use has been made of d3.17b) in the second line. The 
result follows. 



Appendix B: Square of the Dirac distribution 

One possible definition of distribution products, due to 
Colombeau |p70t l7Tll . is as equivalence classes, whose repre- 
sentative chosen for calculations must be inferred from the 
context. This concept allows one to work properly with ob- 
jects such as S 2 , which is proportional to 5, but with non- 
standard (infinite) proportionality constant ll7Tll . 

In the one-dimensional case for instance, take 5(x) — 
linio-^o S a (x) where S a (x) = e~^ x l a ^ I 2 / '(V^cr) is a Gaus- 
sian delta-sequence. Then J dx5 2 (x) = l/(cr-\/47r). This al- 
lows one to set a5 2 (x) = qS(x) with q = 1/^/4tt, as a — » 0. 
A different choice of delta-sequence would lead to some other 
q. This quantity therefore depends on the choice of representa- 
tion of S(x), motivated by the physical nature of the problem 
considered. The d-dimensional generalization of this argu- 
ment leads to d3.27| ). 

We note that, in principle, a more fundamental treatment 
of the d-dimensional case would require acknowledging that 
the Dirac contribution to the Green function (13.2) stems from 
a representation (using a superscript to emphasize the dimen- 
sion) UM 

^)(r) = lim<5f(r), 5^(r) = f^, (Bl) 

where S(r — r/), with r\ identical to that in the principal 
value prescription in Eq. (13.2) . materializes the surface of the 
Lorentz cavity associated to the exclusion volume between in- 
teracting polarizable elements [46]. Doing so would provide 
the q of Eq. (13.27) as a function of r\, but also inevitably intro- 
duce other arbitrary constants, leading to unnecessary compli- 
cations. 
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